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Summary 

Climate changes traditionally have been detected from long series of observations and long after 
they happened. The 8 "inverse sequential" monitoring procedure is designed to detect changes as soon as 
£ey STrequency distribution parameters) are elated both from the most recent existing « of 
observations and from the same set augmented by 1 ,2, • • • j new observations. Indl ^* v ^3P^ 
itv products ("likelihoods") are then calculated which yield probabdiues for erroneously accepting the 
existing parameter(s) as valid for the augmented data set and vice versa. A parameter change is signaled 
when these probabilities (or a more convenient and robust compound "no change probabi 1 ^ s w 
progressive decrease. New parameters are then estimated from the new observations done to restart th 
procedure. The detailed algebra is developed and tested for Gaussian means and variances, Poisson an 
chi-square means, and linear or exponential trends; a comprehensive and interactive Fortran program i 

provided in the appendix. 
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1. Introduction 

The detection of changes in a developing time series requires some idea of what form the changes 
are likely to take. When the nature of the forcing is known, filters can be designed that will show their 
effects most clearly (Kim and North, 1991), but that knowledge is often not available in the geophysical 
sciences. However, many time series are made up of irregular-length sections each of which differs from 
its neighbors in one or more of the parameters that define its signal and noise characteristics. As long as 
its parameters remain unchanged, an individual section can then be said to be in "statistical control" 
(Shewhart, 1939). 

There exists considerable evidence that this concept is realistic in many geophysical contexts. 
Examples are variables exhibiting the "Hurst phenomenon" much discussed in hydrology (e.g., Klemes, 
1974), and atmospheric circulation patterns (Toth, 1992). With its minimum of arbitrary assumptions, the 
concept of statistical control suggests a general monitoring approach that registers the length and end of 
each controlled "regime", together with the new parameter values. The magnitude of changes in geophy- 
sical parameters cannot be anticipated, but their surveillance might use a probability for regarding the 
parameters established from existing observations as significantly changed by the addition of one or more 
new observations. 

Such a "sequential" use of accruing information was pioneered by Wald (1947) and has developed 
into a large special field of statistics (c.f. e.g., Gosh, 1988) which includes a range of procedures utilizing 
cumulative sums ("cusum" techniques; e.g., Goel, 1982). The typical outcome in the simplest situation is 
a decision, with prescribed error probabilities, to accept one of two specified parameter values, or to con- 
tinue sampling. 

The "inverse" sequential approach here presented instead progressively determines "no change" pro- 
babilities for parameter estimates based, respectively, on the accrued data and on the same data aug- 
mented by one or several new observations. A parameter change is then signaled when these probabilities 
begin decreasing to small values. 
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The basic approach is developed in the next section and formulated in Section 3 for the parameters 
of Gaussian, Poisson, and chi-square distributions, as well as for linear and exponential trends. The tech- 
nique is illustrated there with constructed examples, and is described as a systematic and comprehensive 
procedure in Section 4 with reference to an interactive Fortran program, reproduced as code in the 
appendix. 

2. Theory 

Consider a series of m observations *, , i = 1 ,2 • • • m , to which further j observations are added 
(j = 1 ,2, • • ■ ). For a parameter 0 (such as mean, variance, trend, etc.) the first m values yield an optimum 
estimate 9 m which the augmented set of observations changes to 0 m+7 . Writing the corresponding proba- 
bilities of individual x t as p (x , ; 0 ), the likelihood function of n observations is 

L{n\ Q) = p(x l -, 0) p(x 2 , 0) p(x n ;Q). (1) 

With n = m or m + j and 0 = 0 m or Q m+J we have four different likelihoods: 

L l = L{m;Q m )-,L 2 = L(m +j ; 9 m+; ) ; L 3 = L(m +j ; 0 m ) ; L 4 = L(m; B m+j ) . (2) 

Now the form of L shows that the likelihoods decrease systematically with increasing sample size. 
Those for the initial data can be made comparable to those for the augmented data by multiplying the 
L(m; 0) by some factor c(m) and the L(m +j ; 0) by c(m +j). Furthermore the sum of the adjusted 
likelihoods, c(m)L x + c{m + j )L 3 represents the probability that 0 m is valid for either the initial m data 
or the augmented set of m + j . Since these are taken to be the only choices, that probability is one, the 
same applies to the sum c(m + j )L 2 + c(m )L 4 . Denoting c(m + j )Lj by a and c ( m )L$ by P, the other 
adjusted probabilities become c (m )L 1 = 1 - cx and c{m +j )L 2 - 1 — p. 

(In the terminology of the theory of hypotheses (e.g. Hoel, 1966), a represents the "type I error" 
probability for not accepting 0 m , even though true, for the augmented sample. The probability P is that of 
rejecting 0 m+; though true for the initial data; alternatively, it is the "type II" error probability of prefer- 
ring 0 m , though false, for the initial sample). 
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The probabilities a and P are calculated from the factor-free likelihood ratios: 


q(m) = 


ki 

L , 


P 

1 -a ’ 


(3a) 


Solving for a and p yields 


f-2 

q(m+j) = — 
l 3 


1-P 

a 


(1 ~q m ) 

(Qm+j ~ Qm ) 


(3b) 


(4a) 


(?m+y q<n Qm ) 

(<?m +j ~<lm) 


(4b) 


With a definite change of control from 0 m to Q m+J , both probabilities in due course decrease to small 
values. While the existing regime continues, the variability of the likelihood function and rounding errors 
can raise the likelihood ratio q (m ) to unrealistic values larger than unity and similarly reduce q ( m +j ) to 
values below one. To avoid probabilities that are negative or larger than 1, such q values must be 
replaced by 1, implying equality of the likelihoods involved. 

For the inverse sequential monitoring operation, several combined quantities suggest themselves as 
more stable than a and p: 

i) a + P, the probability that either 0(m ) is valid for the sample of m +j or 0(m + j ) for the sample 
of m ; 

ii) a • P, the probability that both these statements are true; 

iii) a compound "no-change" probability y which will be used for the illustrations in Section 3, and is 
defined by writing 

qim ±J) LM = (l-a)d-P) (1-Y) 2 (5) 

q{m) L 3 L 4 aP y 2 


so that 
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Y=(i+V5T‘ • (6) 

The probability y falls between 0 and 0.5 as long as q(m+j) > q(m) and lies between the arithmetic 
and geometric averages of a and p, as can be shown by alternatively substituting these averages for a and 
P in the original form of (5), i.e. in 


[ 1 - (a+p) + ap] 

ap 


When a = p = (a + p)/2, equation (7) becomes 


(7) 


1 - (a + P) + ap + — 

2 + p 2 ‘ 
2 


a a 2 + P 2 

[° P + _ T^J 



(8a) 


This shows that the numerator N and the denominator D of Q both have been increased by 
e = (a 2 + p 2 yi > 0. Since Q> 1 (i.e., N > D), then Q=N/D>Q\= (N + e)/(£> + e) since 
ND + Ne > ND + D e, or N > D, the initial condition. 


Again with a = P = (aP ) , equation (7) becomes 


02 = 


1 - 2(ap y 2 + ap 

ap 


(8b) 


so that Q 2 - Q =- 2(aP) ,/? + (a + P) > 0 since a + P > 2(ap/ /2 ; this can be seen by squaring both sides 
giving 


(a + p) 2 + a 2 + P 2 + 2ap > 4ap , or (a - P) 2 > 0 . 
Finally, with Q 2 > Q > Q \ , 




-1 


~ 7 geometric < Y ~ 




^ 7 arithmetic 


r 'S . 




(9) 


(Equations (3a) and (3b) have the form of the decision limits of Wald’s (1947) "sequential probabil- 
ity ratio test (SPRT)". Log e Q can then be interpreted as the logarithmic width of the indecision region of 
a SPRT in which log e q{m + j) defines the upper decision limit, and log e ^(m) the lower decision limit, 
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respectively). 

A change in parameter(s) lowers y at a rate that increases with the change magnitude but decreases 
with an increase in the number of observations before the change. This is further discussed and illustrated 
in Section 3.1; it suggests using a moving base period, or restarting the procedure with new base values 
after some interval in which y shows no clear descending trend. 

3. Formulae and examples 

The formulae give the basic probability p in the likelihood functions for m and m +j observations, 
and the likelihood ratios q (m ) and q (m +j ) used to calculate the probabilities a, p, and y from equations 
(4) and (6) in Section 2. For simplicity, subscripts will be used to indicate the number of values used for 
parameter estimates, and bracketed symbols for the numbers used to calculate the likelihoods and their 
ratios. Thus L (m ; 0 m ) = L , becomes L m (m), L(m + , 6 m ) = L 3 = L m (m + j ), etc. 

The examples in this section use constructed data with known properties and illustrate how the 
detailed properties of the inverse sequential procedure will be established by more extensive calculations 
using different base lengths and magnitudes of parameter changes. 

3.1 Poisson mean (= variance) 

This case is rather simple because the basic probability 

(3.1.1) 

x \exp (* ) 

has only a single parameter, the mean number * of occurrences. The logarithmic likelihood functions are 

m _ 

log e L m (m ) = mT m log e T m - £ log,* ! - mx m ; 

1 

m _ 

log e L m+ j(rn)=nix m log e x m+J - £ log,* ! -mx m+J ; 


(3.1.2b) 



Table 1 : Detection of change in Poisson mean (= variance). 


The 5 base values are drawn at random from a Poisson population with mean 5. These and 
another 10 values from the same population are used in test III. Test I and test II each use the 
same base values as in test HI and 10 values from Poisson populations with means 3 and 7, 
respectively. Base values used are 6, 1, 6, 6, 4 and have a sample mean of 4.6. 


Test I Test H Test III 


no-change 

probability 

no-change 

probability 

no-change 

probability 

2 

0.47 

5 

0.50 

3 

0.49 

1 

0.36 

6 

0.49 

2 

0.44 

6 

0.43 

9 

0.40 

4 

0.43 

1 

0.33 

4 

0.43 

6 

0.47 

2 

0.25 

9 

0.32 

5 

0.48 

2 

0.18 

1 

0.43 

6 

0.50 

5 

0.22 

6 

0.41 

7 

0.50 

2 

0.16 

10 

0.28 

1 

0.49 

1 

0.09 

5 

0.29 

8 

0.50 

0 

0.04 

13 

0.11 

5 

0.50 
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m+j 

log e L m+J (m +j) = (m+j)x m+J \og e x m+] - X log,x!-(m +j)x m+J ; 

i 


(3.1.2c) 


m+y 

log, L m (m +j) = (m +j )x m +J log, - X ,0 & x\-{m +j )x„ 

i 


(3.1. 2d) 


resulting in the likelihood ratios 

q(m) = exp 


f - , X m +) 

m Um log, _ 

1 

1 

l Xfn 

JJ 


(3.1.3) 


q ( m +j ) = exp 


(m + y)j* m+> log, -^—-(x^ ~~ x m ) I 


(3.1.4) 


so that 


Q = exp 


j (m +j )x m +j - mx m log, - j (*m +j ~ x m ) 


(3.1.5) 


Examples 

Three tests were conducted to demonstrate how the no-change probability is used to detect changes 
in Poisson mean. Table 1 lists three series each containing ten values drawn at random from Poisson 
populations with means of 3, 7, and 5 respectively. All three tests use the same five base values, which 
have a sample mean of 4.6. In tests I and II, progressive decreases occur in the no-change probability 
suggesting a change in the Poisson mean. This is expected since the values in each series were drawn 
from different population means than those of the base period. In test III, the no-change probability does 
not steadily decrease indicating that the series mean is not significantly different from the base period 
mean. This too is expected since the series values and the base values were drawn from the same popula- 
tion. 

The relatively simple form of equation (3.1.5) makes it possible to explore the dependence of the 
no-change probability on change magnitude and base length for an average Poisson sample. Let , 


No-chang* probability No-chang* probability 


Figure 1. No-change probabilities for change in Poisson mean at m=5 or 10. 
Initial mean x; new mean x’ = x + Ax; g = Ax / x 
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i = 1,2 - • • m now represent the average of many samples drawn from a Poisson distribution with mean 
x, so that each x , equals x , and add further j values each equal to x' from a Poisson distribution with 
mean x With Ax — x'-x, Equation (3.1 .5) can be modified to read 


log, Qp = x m 


- 


■ 

f. , Axl , 

Ax 

. Ax 

\j +(m+j) —Hog, 

1 + — 


l xm J 

Xfn j 

x m 


(3.1.6) 


Now the mean x m is simply equal to x , while 


Xm+j 


_ mx + jT _ (m+jytm + jx m (x'-x)/x 
m+j m+j 


(3.1.7) 


With (x'-x )/x = g and Ax/x m = jg /(m+j) = v, say, (3.1.6) becomes 


log, Q P = x„ 


[j +(m +j)vHog e (l + v)-yv 


which together with equation (6) leads to the no-change probability 

l-i 


Y= 


1 + 


(3.1.8) 


(3.1.9) 


as function of x , m , j , and g = Ax/ x m 

Equations (3.1.8) and (3.1.9) have been evaluated for three different x =x m and five different values 
of g . The results are shown in Figure la for three different x' and in Figure lb for a single value x = 5, 
for j — 1 through 10. 

The two broken lines in Figure lb show that a longer base (larger m ) will slow the response of the 
no-change probability to the change in mean, but leaves the curves essentially unaltered in shape. The 
decline in the no-change probability starts at the change of mean and accelerates down to values of the 
order of 0.3— 0.2 before becoming more and more gradual as small probabilities are approached. Beyond 
this general feature, the simple argument here used reveals little about the rate of decline except that it 
changes markedly with both the base mean and the means ratio g . More specific features remain to be 
established by numerical experimentation. 
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3.2 Inverse sequential formulae for Gaussian means and variances 


The basic probability for the Gaussian distribution, 


p = (27 to 2 ) 1/2 exp 


- (x - It) 2 


2a 2 


(3.2.1) 


involves two parameters, p and a 2 , which can only be tested jointly since in the present context neither is 
known a priori. We use the sample mean as an estimate for the population mean p, and 
a 2 = [(n/(n - l)]j 2 , where s 2 is the sample variance and n(=m or m + j) is the number of values used. 
The likelihood products (equation 1) are converted to sums by taking logarithms; thus, 


, , , x m , m , 2 £ (x - p m ) 2 

log, L m (m ) — log, 271 log, G m X _ 2 

1 l 2 O m 


(3.2.2) 


with XU ” P) 2 = (m - l)Om, last term reduces to [~(m - l)/2]. 


Proceeding in the same way for L 4 = L m +7 (m ) leads to 


log, L m+ j{m ) — log, 2 ti log, o m+ y X 2 


m 


" (x-\i m +j ) 2 


2o, 


m+j 


(3.2.3a) 


with \i m+ j - p m = Ap, the numerator of the last term can be written as 


- - (Hm + Ap)) 2 = - “ M-m Y ~ X(- 2xA\i + 2p m Ap + Ap 2 ) 

i l i 


= - (m - 1 )0 2 - 2m p m Ap + 2m p m Ap - m (Ap) 2 , 


so that equation (3.2.3a) becomes 


(jn ) — “ ^ log27i - ^ log^G m 4 .j 2 a 2 


m 


m 


J m+j 


2 Om+j 


(3.2.3b) 


Finally, subtracting equation (3.2.2) from equation (3.2.3b) gives the log likelihood ratio 
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log, q(m) = log. 


L m +j( m ) m m — 1 


L m {m) 2 


= V lo g. 


J m+j 


1 - 


J 


m (A|i) 2 

2Gm+j 


Proceeding the same way for the augmented set of m+j values yields first 


log* L m+J (rn +j) = - log, 2n - log, o 2 + , - £ 


m+j(x -jW;) 2 

2o 2 + , 


1 ^m+j 


m+j 


where the last term, with £ (x - |i m + j ) 2 = (m+j - 1 ) O^+j > reduces to - (m + j - 1 )/2. Next, 


m+j , - m+j 2 "it 2 — (M^m+y Ap)) 2 

log, L m (m +j ) — ^ ^®g? 2tt 2 20 2 


Expanding the last term as before yields 


, r / x M+j , „ m+j , t ( w+ 7 _ 0O/n+; 

log, L m (m +/ ) = — log, 27t — log, o m 

2 2 a~ 


(m+j) Ap 2 

2o 2 


Subtracting equation (3.2.6b) from equation (3.2.5) we obtain the second log likelihood ratio 


, ✓ m+j ’ . m+j - 1 

log, q( m +j) = ^ log, 2 *" — ^ — 

2 °m + , 2 


r„ 2 


■’m+7 


- 1 


(m+jXAp) 2 

2o 2 


The final formulae therefore are 


q(m) = exp 


mlog e 


m -1 
4 - 


J m+j 


1 - 


J 


/n 


2<J 2 +; 


(M^m+y M-m )~ 


and 


q ( m +j ) = exp 


(m +j )log. 


°m+ 1 


m+/-l 


J m+j 


- 1 


^-(p„, +; -p m ) 2 

2a„ 


The first two exponents in each formula reflect solely changes in variance, while the third 


(3.2.4) 

(3.2.5) 

(3.2.6a) 

(3.2.6b) 

(3.2.7) 

(3.2.8a) 

(3.2.8b) 

exponent 


depends primarily on changes in the mean. 



Table 2: Detection of change in Gaussian mean and variance. 


The 5 base values for tests I and II are drawn from a Gaussian population with mean 5 (sample 
mean: 3.64) and variance 25 (sample variance: 71 .74). The 5 base values for test III come from 
a Gaussian population with mean 7.5 (sample mean: 8.02) and variance 6.25 (sample variance: 
4.20). Test I uses 3 values from the first series followed by 7 values from the second. The 10 
values for test II all come from the first series, and those for test HI from the second series. 


Test I 

Test II 

Test in 

- 5.0 

base values 

- 5.0 

base values 

10.2 

base values 

- 1.3 


- 1.3 


9.2 


10.2 


10.2 


6.8 


14.9 


14.9 


8.8 


- 0.6 


- 0.6 


5.1 


2.6 

no-change 

probability 

0.49 

2.6 

no-change 

probability 

0.49 

6.7 

no-change 

probability 

0.49 

6.6 

0.46 

6.6 

0.46 

6.9 

0.45 

2.5 

0.41 

2.5 

0.41 

5.1 

0.40 

10.2 

0.39 

12.4 

0.41 

9.7 

0.46 

9.2 

0.35 

10.5 

0.37 

14.0 

0.35 

6.8 

0.29 

8.9 

0.32 

10.4 

0.33 

8.8 

0.24 

8.7 

0.27 

6.1 

0.36 

5.1 

0.19 

7.7 

0.21 

6.5 

0.39 

6.7 

0.14 

8.6 

0.17 

8.6 

0.42 

6.9 

0.10 

8.8 

0.12 

6.8 

0.44 
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Examples 

Table 2 uses values drawn at random from two Gaussian populations with different means and vari- 
ances (p = 5, a 2 = .25, and p = 7.5, c 2 = 6.25, respectively). Test I uses the first 5 values of the first series 
as base and continues with the next three values of the same series before introducing those of the second 
series. A steady decrease in the no-change probabilities starts with the third value after the change. By 
contrast in test HI, the second series without major change provides a definite no-change signal. 

In test II a regime change is suggested even though the series used was designed without such a 
change. The mean of the initial 5 "base" values (3.64) is considerably smaller than the population mean 
(5.0) while the entire series of 15 values has a mean (6.37) that is considerably larger than the population 
mean. Acting on the sample information alone (all the inverse sequential procedure is designed to do), 
the test therefore quite properly suggested a significant change from the base parameter. 

3 3 Trends (least-square regression) 

Observations made at equally spaced times t = 1,2 ■ • ■ n (n =m or m + j) are represented by 

y = A +Bt +e . (3.3.1) 

This also covers the case of exponential regression when y = log x . The residuals e are assumed to be 
normally distributed with zero mean and variance c 2 . Sample estimates of the regression coefficients A 
and B satisfying least-square requirements are 

XO -y)(t - T) 

a=y+bT,b=— , (3.3.2) 

ia-t ) 2 

i 

where 

7 = (W + 1) ■£«-T) 2 = n(n2 -~ Jl . 

1 2 ’Y ’ 12 

The regression estimate y' for a given i' is then y' = y +b(t' -T), and the corresponding residual 
e t = y', -y' t has a Gaussian distribution with zero mean and variance (see e.g., Anderson and Bancroft, 
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1952, Section 12.2) 


v 2 = 
s u 


n-2 


£ (y« -yt) 
1= 1 



' 

' . in + 1) 

2- 

rt(/t 2 -l) n 


” +1 +- 

2 \ 


12 


n 

n(rt 2 - 1) 




12 



(3.3.3) 


The general form of the likelihood functions defined by the n residuals e, 
t = 1,2, ■ ■ ■ n(-m or m + y ) is 


l * 


L = exp X 

(=i 


O'/ -//) 


25/. 


(3.3.4) 


For the inverse sequential detection of trend changes equation (3.3.4) takes the following forms: 

* O' -/m )/ 2 


log e L m =log e L! = - 


m -2 


SS(m)-6„d(m)j 

*=1 


j* ^ 

m + 1 

t — 



/n + 1 , 

2 

- 


(3.3.5) 


log, L m+J {m+j)- log, L 2 = ■ 


m + j - 2 


2 SS(m + j)-b 2 +J d(m + j) 


m 


m+j 

• I 

/=! 


<f(m) 


O' -fm+j)! 


m+J_ + 1 
m+j 



d(m +j ) 


(3.3.6) 


log,L m (m + ; ) = log, L 3 = - 


m+j - 2 

^ O'-Jm) 2 

2 

SS (m + j)-b£d{m +j) 

. 

z- 

f 1 

m + j + 1 
; 

2 


w+j + l 4 l _ 

m+j dim +j) 


(3.3.7) 
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log* L m+ j (m ) = log* Z.4 = — 


m — 2 


(y -y> 


m+} 


), 2 


SS{m)-b 2 +J d(m)^ 

■ 2a ~ 

1 = 1 

r 

* ' 

m + 1 

2 



m + 1 

2 

t J 




m 

dim) 



; (3.3.8) 


with 


SS(n)=£(y-y ) 2 ; 

/=i 

& =a n + b n t 

Finally as before 

q(m) = e\pilog e L 4 -log e L l ) , 

and (3.3.10) 

q(m +;) = exp(log, L 2 - log, L 3 ) . 

Examples 

Table 3 gives the results of four tests to demonstrate the inverse sequential detection of changes in 
linear trend (regression coefficient b in y = a + bt). The test III series from Table 2 is used. This series 
starts with declining values but settles down to a negligible trend for the entire sample 15 values 
(h =-0.02; variance s 2 = 5.88). Positive and negative trends of b = ±s/ 6 and b =±5/3 are then 
imposed on the last ten y -values in the series. 

The no-change probability starts decreasing in each case as soon as the imposed trend produces a 
distinct change from that of the base series. By contrast, when the last 10 values of the base series are 
tested unchanged all the no-change probabilities (not shown) remain above 0.45. 


d(n) = 


n(n 2 - 1) 


12 

n = m or m + j 


(3.3.9) 


3.4 Chi-square 

The variances 5 2 of samples from a Gaussian population with variance o 2 define a chi-square variate 


x 2 = 


hs 2 
a 2 ’ 


(3.4.1) 



Table 3: Detection of changes in linear trend 


Base series used is that of test III in table 2, with the following sample 
parameters: mean 8.10; variance = 5.88; regression 
coefficients a = 8.27, b = - 0.021 

1 ) trend b = ± 0.4 imposed on the 1 0 y -values from test III, table 2 


y 

positive trend 
b 

(progressive) 

no-change 

probability 

y 

negative trend 

b no-change 

(progressive) probability 

7.1 

-0.04 

0.50 

6.3 

-0.06 

0.49 

7.7 

-0.04 

0.50 

6.1 

-0.08 

0.48 

6.6 

-0.06 

0.49 

4.2 

-0.14 

0.42 

11.3 

0.00 

0.50 

8.1 

-0.11 

0.45 

16.3 

0.12 

0.38 

12.3 

-0.03 

0.50 

12.8 

0.16 

0.29 

8.0 

-0.02 

0.50 

8.9 

0.14 

0.30 

3.3 

-0.07 

0.48 

9.7 

0.13 

0.29 

3.3 

-0.11 

0.44 

12.2 

0.15 

0.23 

5.0 

-0.13 

0.41 

10.8 

0.15 

0.20 

2.8 

-0.15 

0.35 


2) trend b = ± 0.8 imposed on the 10 y -values from test III, table 2 


positive trend negative trend 


y 

b 

(progressive) 

no-change 

probability 

y 

b 

(progressive) 

no-change 

probability 

7.5 

-0.03 

0.50 

5.9 

-0.07 

0.49 

8.5 

-0.02 

0.50 

5.3 

-0.11 

0.45 

7.8 

-0.02 

0.50 

3.0 

-0.17 

0.34 

12.9 

0.06 

0.45 

6.5 

-0.17 

0.35 

18.3 

0.19 

0.23 

10.3 

-0.10 

0.44 

15.2 

0.25 

0.10 

5.6 

-0.11 

0.42 

11.7 

0.25 

0.09 

1.0 

-0.18 

0.30 

12.0 

0.25 

0.08 

0.1 

-0.23 

0.17 

15.8 

0.28 

0.03 

1.4 

-0.26 

0.10 

14.8 

0.30 

0.02 

-1.2 

-0.30 

0.04 
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where ft is the number of values in each sample. If these values are independent of one another, x 2 has 
ft-1 degrees of freedom (d.f.). For "coherent” (autocorrelated) series the d.f. number (which also 
represents the mean of the chi-square distribution as well as one half its variance) is reduced to 

X 2 = v = ft - £(ft ) , (3.4.2) 


(Radok, 1992) where 


£(ft ) = 1 + 


2(ft -1) 


r i + 


2 (ft - 2) 


r 2 + , ■ ■ ■ , + -j- r h 


(3.4.3) 


Here the r, are the autocorrelations of observations i values apart, and ft-e(ft) represents the 
number of independent observations in each section, which equals ft-1 when all autocorrelations are 
zero. 


The basic probability for the chi-square distribution is 



■ 

r > " 

p = 

T n r 

V 

2 

< J ^ 


-i 


(X 2 ) (v/2) -> exp 


x: 

2 


(3.4.4) 


Here v is the number of degrees of freedom which equals the mean as well as half the variance of the dis- 
tribution. Then the logarithmic likelihood functions take the form 


logL m (rn) = 


m 


v w log,2-mlog,r 


y m 

2 



m i m 

Iiog e x 2 -7lx 2 ; 

1 Z 1 


(3.4.5a) 


m 


log,L m+J (rn)- v m+j/ log e 2 tn log e T 


v m+j 


v m +j 


- 1 


m i m 

-I}og e x 2 -^ lx 2 ; 

i z i 


(3.4.5b) 


log e L m+ .(m +;') = - 


m+j 


v m+y log e 2-(m +y')Iog e r 


v m+j 


+ 



m+j ] m+j 

Siog f x 2 -7Ex 2 ; 

i z i 


(3.4.5c) 
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I°g e L m (w +;) = - 


__ ™ ±1 

2 


v m log f 2 - (m +j )log f T 


-1 


m+j i m+j 

2>g j 2 -t!x 2 


The likelihood ratios become 


q(m) = 2 


- 0 2 


(V* — Vm*} ) 


2 

J 

^m+/ 


1 m 

exp ~(v m+y -v m )£logz 2 

z i 


and 


<?("! +j) = 2 


m + 2 v «+y) 



f 


r 

V m 


2 



v m+j 


m+j 


exp 


y(v m+ j - V m )£logx 2 


so that 


logf Q = exp 


2 ^m+/ ) 2 + y iog e 


v m+> -v m m+; 


v m+j 


£iog,r 

/n + 1 


The gamma functions are evaluated with the Euler relation 


where 


nzr'=(z)e« n 

i=i 


i + 


, -z/i 




(3.4.5d) 


(3.4.6) 


(3.4.7) 


(3.4.8) 


(3.4.9) 


c = lim 


= 0.5772156649 - • • 


(3.4.10) 
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Examples 

Tests of the inverse sequential procedure for detecting changes in coherence use three series of 
values drawn from a random independent Gaussian data set. As described by Radok (1992), these data 
were rendered coherent by applying 3-term and 7-term moving averages that produced for different lags X 
autocorrelations of magnitude r^ = (M - X)/M , where M is the length of the moving average. According 
to equations (3.4.2) and (3.4.3), the chi-square distributions of sample variances computed from sets of 5 
successive values should have 2.53 and 1.14 degrees of freedom (d.f.), respectively, compared with 4 for 
the 5-term sample variances from the original series. 

Table 4 shows the chi-square values used for testing the ability of the inverse sequential procedure 
to detect changes in coherence. Test I uses the first 5 values of the series with 2.53 d.f. as base and con- 
tinues with the chi-square values of the non-coherent series; this implies a nominal change from 2.53 to 4 
degrees of freedom. Test II after the same 5 initial values, continues with the chi-square values derived 
from the series of 7-term moving averages, implying a nominal change from 2.53 to 1.14 degrees of free- 
dom. Finally test III uses the chi-square values of the 2.53 d.f. series throughout. 

In test I the loss of coherence (increase in degrees of freedom) is shown by a slow but ultimately 
clear decrease in the no-change probability. In test II the increase in coherence does not lower most of 
the no-change probabilities below 0.40; a larger base period might have rendered the test more sensitive 
for this small signal. The probabilities of test II are similar to those of test III in which only one series of 
chi-square values was used to simulate absence of change. 

4. Implementation 

The inverse sequential procedure is designed to detect changes in statistical control from a few new 
observations added progressively to a representative sample drawn from the most recent controlled data 
regime. As an initialization (which also establishes the presence or absence of such regimes in the exist- 
ing data) the procedure is applied to the full available data set. For this, progressive means, variances. 



Table 4: Detection of a change in chi-square degrees of freedom 
(d.f. = mean = one half variance) 

All tests use as base the following 5 values drawn from a chi-square series constructed with 2.53 
degrees of freedom (for details cf. section 3.4): 

2.72, 3.88, 2.12, 0.18, 1.27 sample mean (d.f.) 2.034 

Another 15 values drawn from the same series are used for test m. Test I uses 15 values from a 
chi-square series constructed to have 4 d.f., while the data for Test II come from a series con- 
structed to have 1.14 d.f. 


Test I Test II Test ffl 



no-change 

probability 


no-change 

probability 


no-change 

probability 

4.31 

0.47 

0.40 

0.48 

0.79 

0.50 

5.33 

0.39 

5.00 

0.50 

4.04 

0.50 

1.07 

0.42 

0.10 

0.49 

1.36 

0.50 

5.03 

0.34 

0.69 

0.47 

4.67 

0.42 

3.21 

0.30 

0.08 

0.40 

2.03 

0.47 

0.57 

0.37 

1.11 

0.39 

2.14 

0.46 

3.08 

0.34 

0.84 

0.38 

3.16 

0.42 

1.54 

0.35 

3.03 

0.42 

0.59 

0.46 

3.19 

0.32 

1.02 

0.42 

2.26 

0.45 

0.35 

0.38 

1.20 

0.43 

2.96 

0.43 

1.70 

0.39 

2.35 

0.45 

3.86 

0.39 

2.35 

0.37 

1.80 

0.47 

1.27 

0.40 

3.85 

0.33 

0.50 

0.45 

1.49 

0.41 

14.1 

0.18 

0.44 

0.44 

8.74 

0.29 

12.3 

0.10 

0.42 

0.42 

3.37 

0.26 
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and both linear and exponential regression coefficients are calculated as basic information. In addition, 
the data are divided into successive small subsamples; their means and variances, apart from being used 
for determining changes in coherence (c.f. Section 3.4) provide indications of the underlying probability 
distribution. For Gaussian data these means and variances are independent of one another; a linear depen- 
dence of sample variance on sample mean can be removed (and hence normality created) by taking the 
square root of each observed value, and a dependence of both mean and variance on the sample size is 
eliminated by a logarithmic transformation (for details see Kendall and Stuart, 1966, Chapter 37). 

The Fortran program in the appendix has been designed interactive to allow for the fact that in prac- 
tice some parameters may not have to be considered. For instance, no stationary mean can exist in the 
presence of a clear linear or exponential growth; for discrete rate events a Poisson rather than a normal 
mean is alone relevant, and only a chi-square variable is involved in the test of variance and coherence 
changes described in Section 3.4. 

In order to establish base values of the parameters that do describe the current data and their varia- 
tion, the inverse sequential procedure is applied backward from the most recent observation to find the 
most recent time interval within which the no-change probability remains high. Questions concerning the 
optimum length of such a "base period", and the efficiency of the inverse sequential procedure in detect- 
ing a given parameter change magnitude, will be addressed in systematic experiments with constructed 
sample series during the remainder of the project, together with real-time tests of some of the GEDEX 
data (Olsen and Wamock, 1992; Schiffer and Unimay ar, 1992). 

When several parameters are tested they will in general not lose control at or even near the same 
time. The successive no-change probabilities of a single parameter, as well as concurrent probabilities for 
several different parameters and/or separate series, can be combined, following Fisher (1941, Section 
21.1) as a "fingerprint" of change in the form of a chi-square variate with 2k degrees of freedom 
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k 

l2kd.f. = X(- 2\og e y k ) 
i 

where k is the number of probabilities thus combined. 

Finally, it must be emphasized that the inverse sequential algorithm is intended for the exploration , 
rather than a confirmation , of parameter changes. Flueck and Brown (1993) have shown that an explora- 
tion can be carried out without the panoply of rigorous statistical procedures needed for a confirmation. 
Even so the full properties of exploratory parameters such as the no-change probability y deserve to be 
clarified with numerical experiments planned for the remainder of this project. 

Acknowledgement: Support for this work has been provided by NASA Grant NAGW-2706. Partial sup- 
port for the second author was provided by NOAA’s Climate and Global Change Program. 
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Appendix 

Program "SEQUITOR" is designed to be an interactive program for analysis of Gaussian 
mean and variance, Poisson mean, chi-square (coherence), and linear (or exponential) trend 
changes in a sequential time series. The user typically will receive FORTRAN source code, 
providing an opportunity to make code changes as desired. For example, in the original code, data 
input is assumed to be free format. However, the user may desire to change this to a specific 
format It may also be desirable for the user to add write statements that exclude headings, such 
that the results can then be easily imported into a graphics package. 

An input/output flowchart is included in this appendix. Each square box represents an 
input step by the user, and an oval represents results output. A brief description of the input steps 


follows: 


Enter input filename: This is the input data filename up to 80 characters. 

Enter descriptive title: This a descriptive header of the data and/or the analysis up to 80 characters. 

Enter number of values in series: This is the total number of rows in the input data file. It is 
assumed that the input data file contains a column of x- values (column 1) which represent an index 
or year for example, followed by n columns of y-values containing the actual series for analysis. 

Enter missing value: This program allows for missing values. Enter a unique number (e.g., -999.) 
to represent missing values. 

Enter column number: Input data files may contain multiple y-value columns. This entry should be 
from 1 to n depending upon which y-value column is desired for analysis. The very first column 
in the input data file is considered column 0. 

Enter 0=continue, l=reverse data input order: Often it is desirable to do the sequential monitoring 
analysis beginning with the most recent value and working backwards. This helps identify 
"regimes" in the time series. Enter either a 0 or 1. 


Enter beginning and ending x-range values: Enter values separated by a comma. This range 
corresponds to the x-values in the very first column of the input data file. Since the x-values might 
represent an index or year, examples would be 10,18 or 1985,1992. Note that these values can 
represent a sub-set of the input data file. 


Enter window size for sub-samples: In determining a regime, it is useful to examine smaller sub- 
sets of values. A typical sub-sample might contain 5 values. If the total number of cases in the 
series is not evenly divisible by the window size, the remaining values will be ignored in only the 
sub-sample analysis. 


Enter analysis type: Here there are several options. Entering 1 through 4 places the user in the 
desired sequential analysis routine. Other options include changing the sub-sample size, changing 



the column number, changing the data range, reversing the data order, or simply quitting the 
program. 

Enter number of base period values: Within each analysis routine, the user is prompted for the 
number of base period values. These should typically be small, say 5 to 15 or so. Base period 
results is ouput at this point. 

Enter 0=continue, l=change number of base values: Upon examining the base period results, the 
user is given the option to change the base period size, or continue with the final analysis. 

This program was written interactively because it is intended to be exploratory in nature. 

An attempt was made to allow the user to make changes during the analysis, instead of having to 
restart the program several times. Results are output to the screen and to a file named 
"sequitor.out", which is replaced each time the program is run. 

The program contains minimal comments, but variables are defined at the beginning of each 
subroutine to help the user understand the program. It is intended to use the program in 
conjunction with this progress report, and some attempt at consistency of variable names in relation 
to the formulae has been made. However, it is possible that updates or changes to program will 
occur after this initial release. Hence, users may want to contact the authors for additional 
information. 

Output from a sample analysis is included in this appendix. It is test I from Section 3.2 in 
the progress report. The associated input data are also included. 



Flowchart of input/output for program SEQUITOR 




Enter missing value: 


Enter column number: 



Enter analysis type 

1= Gauss, 2=Pois8on, 3=chi-square, 4= I inear, 
0=change sub-sample size, 6=change column number, 
7=change data range, 8=reverse data order, 9=quit: 


Gauss, Poisson, chi-square, 
or linear trend output 
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1978. 

-5.0 -5.0 10.2 

1979. 

-1.3 -1.3 9.2 

1980. 

10.2 10.2 6.8 

1981. 

14.9 14.9 8.8 

1982. 

-0.6 -0.6 5.1 

1983. 

2.6 2.6 6.7 

1984. 

6.6 6.6 6.9 

1985. 

2.5 2.5 5.1 

1986. 

20.2 12.4 9.7 

1987. 

9.2 10.5 14.0 

1988. 

6.8 8.9 10.4 

1989. 

8.8 8.7 6.1 

1990. 

5.1 7.7 6.5 

1991. 

6.7 8.6 8.6 

1992 . 

6.9 8.8 6.8 



ENTER INPUT FILE NAME: 
gauss . dat 

ENTER DESCRIPTIVE TITLE: 

Detection of change in Gaussian mean and variance. Test I 
INVERSE SEQUENTIAL MONITORING (PROGRAM <SEQUITOR>) 

ENTER NUMBER OF VALUES IN SERIES: 

15 

ENTER MISSING VALUE: 

-999 

ENTER COLUMN NUMBER: 

1 

ENTER 0=CONTINUE, 1=REVERSE DATA INPUT ORDER: 

0 

ENTER BEGINNING AND ENDING X-RANGE VALUES: 

1978,1992 

FULL SAMPLE UNIVARIATE STATISTICS: 

X-VALUE RANGE - 1978.-1992. 

NUMBER OF VALUES = 15 

NUMBER OF MISSING VALUES = 0 

SAMPLE MEAN = 6.240 

SAMPLE VARIANCE = 40.034 

SAMPLE SLOPE = 0.466 


ENTER WINDOW SIZE FOR SUB- SAMPLES: 
5 


JUB-SAMPLE PARAMETERS: 






INDEX 

X-VALUE 

NUMBER OF 

MISSING 




RANGE 

RANGE 

VALUES 

VALUES 

MEAN 

VARIANCE 

CHI-SQUARE 

1- 5 

1978.-1982. 

5 

0 

3.640 

71.713 

7.165 

6- 10 

1983.-1987. 

5 

0 

8.220 

52.852 

5.281 

11- 15 

1988.-1992. 

5 

0 

6.860 

1.723 

0.172 


ENTER ANALYSIS TYPE 

1=GAUSS, 2= POISSON, 3=CHI-SQUARE, 4=LINEAR, 
0=CHANGE SUB- SAMPLE SIZE, 6=CHANGE COLUMN NUMBER, 

7 =CHANGE DATA RANGE, 8=REVERSE DATA ORDER, 9=QUIT: 

1 

+ + 

I TEST FOR CHANGE IN GAUSSIAN MEAN AND VARIANCE I 
+ + 

ENTER NUMBER OF VALUES IN BASE PERIOD: 

5 

BASE PERIOD PARAMETERS: 

X-VALUE RANGE = 1978.-1982. 

NUMBER OF VALUES = 5 

NUMBER OF MISSING VALUES = 0 

BASE MEAN = 3.640 

BASE VARIANCE = 71.713 


ENTER 0=CONTINUE, 1=CHANGE NUMBER 

0 

OF BASE VALUES: 



PROGRESSIVE PARAMETERS: 

X Y 

INDEX OBSERVATIONS MEAN 

VARIANCE 

GAMMA 

DELTA 

GAMMA 
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6 

1983.000 

2.60 

3.467 

57.551 

0.486 


7 

1984.000 

6.60 

3.914 

49.361 

0.456 

0.030 

8 

1985.000 

2.50 

3.737 

42.560 

0.410 

0.047 

9 

1986.000 

20.20 

5.567 

67.353 

0.453 

-0.044 

10 

1987.000 

9.20 

5.930 

61.189 

0.422 

0.031 

11 

1988.000 

6.80 

6.009 

55.139 

0.394 

0.028 

12 

1989.000 

8.80 

6.242 

50.775 

0.351 

0.043 

13 

1990.000 

5.10 

6.154 

46.644 

0.322 

0.029 

14 

1991.000 

6.70 

6.193 

43.078 

0.280 

0.042 

15 

1992.000 

6.90 

6.240 

40.034 

0.236 

0.044 

ENTER 

0 

II 

1 

1=CHANGE 

NUMBER 

OF EASE VALUES: 



ENTER ANALYSIS TYPE 

1=GAUSS, 2= POISSON, 3=CHI-SQUARE, 4=LINEAR, 
0=CHANGE SUB-SAMPLE SIZE, 6=CHANGE COLUMN NUMBER, 
7=CHANGE DATA RANGE, 8=REVERSE DATA ORDER, 9=QUIT: 
9 


END OF SEQUITOR RUN 
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PROGRAM SEQU I TOR 


C AUTHOR: TIMOTHY J. BROWN | 

C | 

C INVERSE SEQUENTIAL PROGRAM | 

C | 

C REVISION HISTORY | 

C LEVEL AUTHOR DATE DESCRIPTION | 

C . 01A . TJB 93/04/28. ORIGINAL VERSION. | 

C I 

c ] 

c 

C SIX SUBROUTINES ARE ATTACHED TO THE MAIN PROGRAM: 

C ' SGAUSS' COMPUTES CHANGE IN GAUSSIAN MEAN AND VARIANCE. 

C 'SPOISS' COMPUTES CHANGE IN POISSION MEAN. 

C ' SCHI' COMPUTES CHANGE IN CHI-SQUARE DEGREES OF FREEDOM. 

C 'SLINEAR' COMPUTES CHANGE IN LINEAR TREND. 

C 'UNIVAR' COMPUTES UNIVARIATE STATISTICS MEAN, VARIANCE, AND SUM. 

C ' RCOEFF ' COMPUTES LLS REGRESSION BO AND B1 COEFFICIENTS. 

C 

C INPUT IS ASSUMED TO BE FREE-FORMAT, BUT USER CAN CHANGE AS DESIRED. 

C 

C THE PARAMETER STATEMENT AND COMMON BLOCK IS LOCATED IN ALL SUBROUTINES. 
C THE USER SHOULD CHANGE ' NDIM' AS REQUIRED. 

C 

C THE FOLLOWING ARRAYS AND VARIABLES ARE USED IN THE COMMON BLOCK: 

C ' FCHI2' CHI-SQUARE VALUE FOR EACH SUB-SAMPLE. 

C ' F INDEX' INDEX VALUE (1, 2, . . .N) ASSOCIATED WITH EACH Y-VALUE . 

C ' FXVAL' INPUT X-VALUES . 

C ' FYVAL' INPUT Y-VALUE S . 

C ' XDATA' WORK ARRAY FOR X-VALUES. 

C ' YDATA' WORK ARRAY FOR Y-VALUES. 

C 

C ' FMISS' NUMBER REPRESENTING MISSING VALUES. 

C 'NDIM' DIMENSION SIZE FOR DATA AND WORK ARRAYS. 

C ' NCASE' NUMBER OF FULL SAMPLE VALUES WITHIN INDEX RANGE. 

C 'OUNIT' OUTPUT UNIT NUMBER. 

C 

C THE FOLLOWING ARRAYS AND VARIABLES ARE USED IN THE MAIN PROGRAM: 

C ' FDATA' HOLDS THE Y-VALUES WHEN THEY ARE INPUT; SHOULD BE 
C DIMENSIONED >= NUMBER OF COLUMNS IN INPUT FILE. 

C ' FXDATA' HOLDS THE ORIGINAL X-VALUES OR REVERSED ORDER VALUES. 

C ' F YDATA' HOLDS THE ORIGINAL Y-VALUES OR REVERSED ORDER VALUES. 

C ' XWORK' WORK ARRAY FOR X-VALUES. 

C ' YWORK' WORK ARRAY FOR Y-VALUES. 

C 

C 'BO' INTERCEPT FROM LLS REGRESSION. 

C 'Bl' SLOPE FROM LLS REGRESSION. 

C ' CFILE' INPUT DATA FILE NAME. 

C ' CTITLE' DESCRIPTIVE TITLE. 

C 'H' NUMBER OF VALUES WITHIN EACH SUB- SAMPLE. 

C ' I ' DO LOOP COUNTER . 

C 'II' INDEX COUNTER. 

C ' IBEG' INDEX COUNTER. 

C ' TEND' INDEX COUNTER. 
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C ' IDIR' 

C 'I TYPE' 

C ' I UNIT' 

C ' J' 

C 'K' 

C 'N' 

C ' NBEG' 

C 'NEND' 

C ' NN f 
C ' NCOL' 

C 
C 

C ' NMISS' 

C ' NPOP' 

C ' NVALS' 

C ' POP VAR' 

C ' S2' 

C ' XBEG' 

C 'XEND' 

C ' YBAR' 

C ' YSUM' 

C ' YVAR' 

C 

PARAMETER (NDIM=150) 

COMMON /WORK/ FXVAL (NDIM) , FYVAL (NDIM) , XDATA(NDIM), 

+ YDATA (NDIM) , FINDEX (NDIM) , FCHI2 (NDIM) , NCASE, 

4- OUNIT, FMISS 

REAL FDATA ( 15 ) 

REAL FXDATA (NDIM) , F YDATA (NDIM) , XWORK (NDIM) , YWORK (NDIM) 
INTEGER H, OUNIT 
CHARACTER* 80 CTITLE, CFILE 

DATA I UNIT, OUNIT / 1, 2 / 

C 

C 

C THIS SECTION REQUESTS THE INPUT INFORMATION, OPENS FILES, INPUTS 
C THE DATA, AND COMPUTES FULL SAMPLE UNIVARIATE STATISTICS. 

C 

WRITE (*, 801) 

READ (*,101) CFILE 

OPEN (IUNIT, FILE=CFILE, STATUS='OLD' ) 

OPEN (OUNIT, FILE=' sequitor . out ' ) 

WRITE ( * , 802) 

READ (*,101) CTITLE 

WRITE (OUNIT, 900) 

WRITE (*, 900) 

WRITE (OUNIT, 901) CTITLE 

WRITE(*, 803) 

READ ( * , * ) NPOP 
WRITE ( * , 804) 

READ ( * , * ) FMISS 


DATA DIRECTION FLAG ( 1=REVERSE DATA ORDER, 0=CONTINUE) . 
ANALYSIS TYPE. 

INPUT UNIT NUMBER. 

DO LOOP COUNTER. 

COUNTER. 

DO LOOP COUNTER. 

BEGINNING INDEX NUMBER FOR INDEX RANGE. 

ENDING INDEX RANGE FOR INDEX RANGE. 

COUNTER. 

COLUMN NUMBER OF Y-VALUES TO BE ANALYZED. 

THIS IS USEFUL FOR FILES CONTAINING MULTIPLE COLUMNS OF DATA. 
X-VALUES ARE ASSUMED TO BE IN COLUMN ONE. 

NUMBER OF MISSING VALUES. 

NUMBER OF POPULATION VALUES. 

NUMBER OF NON-MISSING VALUES. 

POPULATION VARIANCE FROM FULL SAMPLE. 

SUM OF SQUARES IN CHI-SQUARE CALCULATION. 

BEGINNING VALUE OF X-RANGE . 

ENDING VALUE OF X-RANGE. 

MEAN OF Y-VALUES. 

SUM OF Y-VALUES. 

VARIANCE OF Y-VALUES. 


oooon 99 o 9 9 oono 
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3 CONTINUE 

WRITE (*, 805) 

READ ( * , * ) NCOL 
REWIND I UNIT 

INPUT THE DATA AND FILL WORK ARRAYS. REVERSE DATA ORDER IF REQUESTED. 


DO 13 I = 1, NPOP 

READ (IUNIT, *) FXDATA ( I ) , (FDATA ( J) , J=1 , NCOL) 

F INDEX ( I ) = FLOAT (I) 

FYDATA ( I ) = FDATA (NCOL) 

CONTINUE 

CONTINUE 

WRITE ( * , 806) 

READ ( * , * ) IDIR 

IF ( IDIR .EQ. 1 ) THEN 

K = 0 

DO 14 I = NPOP, 1, -1 
K = K + 1 

XWORK(K) = FXDATA (I) 

YWORK(K) = FYDATA (I) 

CONTINUE 

DO 15 I = 1, NPOP 

FXDATA (I) = XWORK(I) 

FYDATA (I) = YWORK(I) 

CONTINUE 

END IF 


FILL WORK ARRAYS WITH DATA WITHIN SELECTED INDEX RANGE AND COLUMN. 

4 CONTINUE 

WRITE (*, 807) 

READ ( * , * ) XBEG, XEND 


DO 18 I = 1, NPOP 


IF ( FXDATA (I) .EQ. XBEG ) NBEG = I 



o o o o o non on on 


Apr 28 16:35 1993 sequitor.f Page 4 


IF ( FXDATA (I) .EQ. XEND ) NEND = I 
18 CONTINUE 

IF ( NBEG .LT. 1 ) THEN 

WRITE (*, 810) 

GOTO 4 

END IF 


IF ( NEND .GT. NPOP ) THEN 

WRITE (*, 810) 

GOTO 4 


END IF 


NCASE = 0 


DO 16 I = 1, NPOP 


IF ( I .GE. NBEG .AND. I . LE . NEND ) THEN 


NCASE = NCASE + 1 


FXVAL (NCASE) 
FYVAL (NCASE) 
XDATA (NCASE) 
YDATA (NCASE) 


= FXDATA (I) 

= FYDATA ( I ) 

= F INDEX (I) 

= FYVAL (NCASE) 


END IF 


16 CONTINUE 


COMPUTE FULL SAMPLE STATISTICS AND OUTPUT RESULTS 

CALL UNIVAR ( NCASE, NVALS, NMISS, YBAR, YVAR, YSUM ) 

CALL RCOEFF ( NCASE, BO, B1 ) 

IRANGE = NEND - NBEG + 1 

WRITE (OUNIT, 902) FXVAL (1), FXVAL (IRANGE) , NVALS, NMISS, YBAR, 
+ YVAR, B1 

WRITE ( * , 902 ) FXVAL ( 1 ) , FXVAL (IRANGE) , NVALS, NMISS, YBAR, 

+ YVAR, B1 

1 CONTINUE 

WRITE (*, 808) 

READ ( * , * ) H 


WRITE (OUNIT, 903) 
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C 

C 

C 


C- 


C- 


C 

c 


c 

c 

c 

c- 


WRITE (*, 903) 

COMPUTE SUB- SAMPLE STATISTICS. 

POPVAR = YVAR 
NN = 0 
K = 0 
N = 0 

IBEG = -(H) + 1 


DO 17 I = NBEG, NEND 

K = K + 1 
IBEG = IBEG + 1 
N = N + 1 


IF ( FYVAL(N) .GT. FMISS ) THEN 

NN = NN + 1 

YD AT A (NN) = FYVAL(N) 

END IF 


IF( K .EQ. H ) THEN 

CALL UN I VAR ( NN, NVALS, NMISS, 
COMPUTE CHI-SQUARE VALUES. 


YBAR, 


YVAR, 


IF ( NVALS .GT. 0 ) THEN 
S2 = 0. 


DO 19 M = 1 , NN 

IF ( YDATA(M) .NE. FMISS ) 

+ S2 - S2 + (YDATA(M) - YBAR) * *2 

19 CONTINUE 


FCHI2 (N) = S2 / POPVAR 
ELSE 

FCHI2 (N) = FMISS 
END IF 

C 

c 

C OUTPUT SUB-SAMPLE STATISTICS. 

C 

II = (I-H) + 1 
I END = IBEG + H - 1 


YSUM ) 
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WRITE (OUNIT, 904) II, I, FXVAL(IBEG), FXVAL(IEND), NVALS, 

+ NMISS, YBAR, YVAR, FCHI2 (N) 

WRITE (*, 904) II, I, FXVAL(IBEG), FXVAL(IEND), NVALS, NMISS, 
+ YBAR, YVAR, FCHI2 (N) 

NN = 0 
K = 0 

END IF 

C 

17 CONTINUE 

C 

C 

C BRANCH OFF TO APPROPRIATE SUBROUTINE, CHANGE SUB-SAMPLE SIZE, 

C CHANGE COLUMN NUMBER, OR STOP PROGRAM. 

C IF CHOOSING CHI-SQUARE, THE NUMBER OF CASES BECOMES THE NUMBER 
C OF SUB-SAMPLE INTERVALS. 

C 

2 CONTINUE 

WRITE <*, 809) 

READ ( * , * ) ITYPE 
C 

IF ( ITYPE .EQ. 0 ) THEN 
GOTO 1 

ELSE IF ( ITYPE .EQ. 1 ) THEN 

WRITE (OUNIT, 1001) 

WRITE (*, 1001) 

CALL SGAUSS 

ELSE IF ( ITYPE .EQ. 2 ) THEN 

WRITE (OUNIT, 1002) 

WRITE (*, 1002) 

CALL SPOISS 

ELSE IF ( ITYPE .EQ. 3 ) THEN 

WRITE (OUNIT, 1003) 

WRITE (*, 1003) 

NCASE = INTERVL 

CALL SCHI 

ELSE IF ( ITYPE .EQ. 4 ) THEN 

WRITE (OUNIT, 1004) 

WRITE (*, 1004) 


CALL SLINEAR 
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ELSE IF ( I TYPE . EQ . 6 ) THEN 
GOTO 3 

ELSE IF ( ITYPE . EQ . 7 ) THEN 
GOTO 4 

ELSE IF ( ITYPE . EQ . 8 ) THEN 
GOTO 5 

ELSE IF ( ITYPE . EQ . 9 ) THEN 

WRITE (OUNIT, 907) 

WRITE (*, 907) 

GOTO 999 

ELSE 

WRITE (*, 905) 

GOTO 2 

END IF 


GOTO 2 

101 FORMAT (A) 

801 FORMAT ( ' ENTER INPUT FILE NAME:') 

802 FORMAT (' ENTER DESCRIPTIVE TITLE:') 

803 FORMAT (' ENTER NUMBER OF VALUES IN SERIES:') 

804 FORMAT ( ' ENTER MISSING VALUE:') 

805 FORMAT (' ENTER COLUMN NUMBER:') 

806 FORMAT (' ENTER 0=CONTINUE, 1=RE VERSE DATA INPUT ORDER:') 

807 FORMAT (' ENTER BEGINNING AND ENDING X-RANGE VALUES:') 

808 FORMAT ( / , ' ENTER WINDOW SIZE FOR SUB-SAMPLES:') 

809 FORMAT ( /' ENTER ANALYSIS TYPE',/, 

+ ' 1=GAUSS, 2=POISSON, 3=CHI-SQUARE, 4=LINEAR, ' , /, 

+ ' 0=CHANGE SUB-SAMPLE SIZE, 6=CHANGE COLUMN NUMBER,',/, 

+ ' 7=CHANGE DATA RANGE, 8=RE VERSE DATA ORDER, 9=QUIT:') 

810 FORMAT (/'RANGE EXCEEDS TOTAL NUMBER OF CASES') 

900 FORMAT (' INVERSE SEQUENTIAL MONITORING (PROGRAM <SEQUITOR> ) ' ) 

901 FORMAT (//, A) 

902 FORMAT (/, 'FULL SAMPLE UNIVARIATE STATISTICS:',/, 

+ 'X-VALUE RANGE = ' , F5 . 0, ' , F5 . 0, /, 

+ 'NUMBER OF VALUES = ',13,/, 

+ 'NUMBER OF MISSING VALUES = ',13,/, 

+ 'SAMPLE MEAN = ',F8.3,/, 

+ 'SAMPLE VARIANCE = ',F8.3,/, 

+ 'SAMPLE SLOPE = ',F8.3) 
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903 FORMAT (/,' SUB- SAMPLE PARAMETERS:',/, 

+ ' INDEX X- VALUE NUMBER OF MISSING' , /, 

+ ' RANGE RANGE VALUES VALUES MEAN VARIANCE' , 

+ ' CHI-SQUARE',/, 

+ ' ==========' ) 


904 FORMAT (13, ' -' , 13, 2X, F5 . 0, ' -' , F5 . 0, 2X, I 9, 2X, 17, 2 (2X,F8.3) , 

+ 2X, F10 . 3) 

905 FORMAT (/,' >>> NOT A VALID SELECTION <«',/) 

907 FORMAT (//' END OF SEQUITOR RUN') 


1001 FORMAT (/,'+ +' , 

+ /,' I TEST FOR CHANGE IN GAUSSIAN MEAN AND VARIANCE I', 

+ /,'+ +' ) 

1002 FORMAT (/,'+ +', 

+ /,' I TEST FOR CHANGE IN POISSON MEAN |', 

+ /,'+ +') 


1003 FORMAT ( 

+ /,'+ + ', 

+ /,' | TEST FOR CHANGE IN CHI-SQUARE DEGREES OF FREEDOM |', 

+ /,'+ +') 


1004 FORMAT (/,'+ +' , 

+ /,' I TEST FOR CHANGE IN LINEAR TREND I', 

+ /,'+ +') 


999 STOP 
END 

csssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss 

SUBROUTINE SGAUSS 


c 

C 

c 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE TO COMPUTE CHANGE IN GAUSSIAN MEAN AND VARIANCE. 


VARIABLES USED IN SUBROUTINE NOT DESCRIBED IN MAIN PROGRAM: 


'BARM' 

' BARMJ' 

' DGAMMA' 
'GAMMA' 
'GAMOLD' 
' J' 

' JJ' 

' Jl' 

' J2' 

'M' 

'MBASE' 

'Q' 

'QM' 

' QMJ' 

' SDM' 

' SDM2' 


MEAN OF BASE PERIOD VALUES. 

PROGRESSIVE MEAN. 

CHANGE IN GAMMA FROM PREVIOUS VALUE. 

GAMMA VALUE. 

PREVIOUS VALUE OF GAMMA. 

PROGRESSIVE VALUE INDEX. 

DO LOOP COUNTER. 

STARTING VALUE FOR DO LOOP COMPUTING PROGRESSIVE VALUES. 
ENDING VALUE FOR DO LOOP COMPUTING PROGRESSIVE VALUES. 
REAL VALUE OF NUMBER OF BASE VALUES. 

NUMBER OF BASE PERIOD VALUES. 

RATIO OF QMJ/QM. 

BASE PERIOD LIKELIHOOD RATIO. 

PROGRESSIVE PERIOD LIKELIHOOD RATIO. 

STANDARD DEVIATION OF BASE PERIOD VARIANCE. 

BASE PERIOD VARIANCE. 


non 
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C 

C 

C 

C 

C 

C 


C 

C 

C 

C 

1 


' SDMJ' PROGRESSIVE STANDARD DEVIATION. 

' SDMJ2 ' PROGRESSIVE VARIANCE. 

'Tl' WORK VARIABLE; FIRST TERM IN EITHER Q (M) OR Q (M+J) EQUATIONS. 

' T2 ' WORK VARIABLE; SECOND TERM IN EITHER Q (M) OR Q (M+ J) EQUATIONS 

' T3 ' WORK VARIABLE; THIRD TERM IN EITHER Q (M) OR Q (M+J) EQUATIONS. 

PARAMETER (NDIM=150) 

COMMON /WORK/ FXVAL (NDIM) , FYVAL (NDIM) , XDATA (NDIM) , 

+ YDATA (NDIM) , F INDEX (NDIM) , FCHI2 (NDIM) , NCASE, 

+ OUNIT, FMISS 

REAL J, M 
INTEGER OUNIT 

FILL WORK ARRAYS WITH BASE PERIOD VALUES, COMPUTE BASE PERIOD 
STATISTICS, AND OUTPUT RESULTS. 

CONTINUE 


WRITE (*, 801) 

READ ( * , * ) MBASE 
C 

DO 20 I = 1, MBASE 

XDATA (I) - F INDEX (I) 

YDATA (I) = FYVAL (I) 

20 CONTINUE 

C 

CALL UNIVARf MBASE, NVALS, NMISS, YBAR, YVAR, YSUM ) 

C 

IF ( NVALS .GT. 0 ) THEN 

WRITE (OUNIT, 1001) FXVAL (1), FXVAL (MBASE ) , NVALS, NMISS, YBAR, 
+ YVAR 

WRITE (*, 1001) FXVAL (1), FXVAL (MBASE ) , NVALS, NMISS, YBAR, 

+ YVAR 

ELSE 

WRITE (*, 1000) 

RETURN 

END IF 

C 

WRITE (*, 802) 

READ ( * , * ) I TYPE 
IF ( ITYPE .EQ. 1 ) GOTO 1 

WRITE (OUNIT, 1002) 

WRITE ( * , 1002) 

INITIALIZE AND COMPUTE FIXED VALUES; BEGIN PROGRESSIVE VALUE LOOP. 

BARM = YBAR 
SDM = SQRT(YVAR) 

SDM2 = YVAR 
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M = FLOAT (MBASE) 

GAMOLD = FMISS 
J1 = MBASE + 1 
J2 = NCASE 
C 

C FILL PROGRESSIVE WORK ARRAYS, COMPUTE PROGRESSIVE STATISTICS, AND 
C OUTPUT RESULTS. 

C 

C 

DO 22 JJ = Jl, J2 


C 

C 

C 

C 


XDATA(JJ) = FINDEX(JJ) 

YDATA(JJ) = FYVAL(JJ) 

CALL UNIVAR ( JJ, NVALS, NMISS, YBAR, YVAR, 
COMPUTE GAMMA. SEE TEXT FOR EQUATION DETAILS. 

IF ( NVALS .GT. 0 ) THEN 


YSUM ) 


BARMJ = YBAR 
SDMJ = SQRT(YVAR) 

SDMJ2 = YVAR 
J = FLOAT (JJ) - M 

T1 = M * ALOG (SDM / SDMJ ) 

T2 = ( (M - 1.) / 2.) * (1. - (SDM2 / SDMJ2)) 

T3 = (M / (2. * SDM2) ) * (BARMJ - BARM) * *2 

QM = EXP (T1 + T2 - T3) 

T1 = (M + J) * ALOG (SDM / SDMJ ) 

T2= ( (M + J — 1 . ) / 2.) * ( (SDMJ2 / SDM2) - 1.) 
T3 = ( (M + J) / (2. * SDM2 ) ) * (BARMJ - BARM)**2 

QMJ = EXP (T1 + T2 + T3) 

Q = QMJ / QM 

GAMMA = 1. / (1. + SQRT(Q)) 

D GAMMA = GAMOLD - GAMMA 


IF ( GAMOLD .EQ. FMISS ) THEN 

WRITE (OUNIT, 1003) JJ, FXVAL(JJ), YDATA(JJ), YBAR, YVAR, 
+ GAMMA, CFLAG 

WRITE (*, 1003) JJ, FXVAL(JJ), YDATA(JJ), YBAR, YVAR, 

+ GAMMA, CFLAG 

ELSE 

WRITE (OUNIT, 1003) JJ, FXVAL(JJ), YDATA(JJ), YBAR, YVAR, 
+ GAMMA, CFLAG, D GAMMA 

WRITE (*, 1003) JJ, FXVAL(JJ), YDATA(JJ), YBAR, YVAR, 

+ GAMMA, CFLAG, DGAMMA 


C 


END IF 
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GAMOLD = GAMMA 
ELSE 

WRITE (OUNIT, 1004) JJ, FXVAL(JJ), YDATA(JJ) 

WRITE (*, 1004) JJ, FXVAL(JJ), YDATA(JJ) 

GAMOLD = FMISS 

END IF 

C 

22 CONTINUE 

C 

WRITE (*, 802) 

READ ( * , * ) ITYPE 

IF ( ITYPE .EQ. 1 ) GOTO 1 

801 FORMAT (/, ' ENTER NUMBER OF VALUES IN BASE PERIOD:') 

802 FORMAT (/,' ENTER 0=CONTINUE, 1=CHANGE NUMBER OF BASE VALUES:') 

1000 FORMAT (' ALL VALUES IN BASE PERIOD MISSING') 

1001 FORMAT (/, 'BASE PERIOD PARAMETERS:',/, 

+ 'X-VALUE RANGE = ' , F5 . 0, ' -' , F5 . 0, /, 

+ 'NUMBER OF VALUES = ',13,/, 

+ 'NUMBER OF MISSING VALUES = ',13,/, 

+ 'BASE MEAN = ',F8.3,/, 

+ 'BASE VARIANCE = ',F8.3,/) 

1002 FORMAT (/, 'PROGRESSIVE PARAMETERS:',/, 

+ ' X Y 

+ ' INDEX OBSERVATIONS MEAN VARIANCE 

+ ' ===== ================ ======== ======== 

1003 FORMAT (15, 2X, F8 . 3, IX, F7 .2, 2 (2X, F8 . 3) , 2X,F6. 3, 2X,F6.3) 

1004 FORMAT (15, 2X, F8 . 3, IX, F7 .2) 

RETURN 

END 

C 

SUBROUTINE SPOISS 
C 

C SUBROUTINE TO COMPUTE CHANGE IN POISSON MEAN. 

C 

C VARIABLES USED IN SUBROUTINE NOT DESCRIBED IN MAIN PROGRAM: 

C 'BARM' MEAN OF BASE PERIOD VALUES. 

C ' BARMJ' PROGRESSIVE MEAN. 

C ' DGAMMA' CHANGE IN GAMMA FROM PREVIOUS VALUE. 

C 'GAMMA' GAMMA VALUE. 

C 'GAMOLD' PREVIOUS VALUE OF GAMMA. 

C 'J' PROGRESSIVE VALUE INDEX. 

C 'JJ' DO LOOP COUNTER. 

C 'Jl' STARTING VALUE FOR DO LOOP COMPUTING PROGRESSIVE VALUES. 

C ' J2' ENDING VALUE FOR DO LOOP COMPUTING PROGRESSIVE VALUES. 

C 'M' REAL VALUE OF NUMBER OF BASE VALUES. 

C 'MBASE' NUMBER OF BASE PERIOD VALUES. 


DELTA' , /, 
GAMMA GAMMA',/, 



ooo o on n noon 
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C 'Q' COMBINED LIKELIHOOD RATIOS. 

C 'Tl' WORK VARIABLE; FIRST TERM IN EITHER Q EQUATION. 

C ' T2' WORK VARIABLE; SECOND TERM IN EITHER Q EQUATION. 

C ' T3' WORK VARIABLE; THIRD TERM IN EITHER Q EQUATION. 

C 

PARAMETER (NDIM=150) 

COMMON /WORK/ FXVAL (NDIM) , FYVAL (NDIM) , XDATA (NDIM) , 

+ YDATA (NDIM) , FINDEX (NDIM) , FCHI2 (NDIM) , NCASE, 

+ OUNIT, FMISS 

REAL J, M 
INTEGER OUNIT 

FILL WORK ARRAYS WITH BASE PERIOD VALUES, COMPUTE BASE PERIOD 
STATISTICS, AND OUTPUT RESULTS. 

1 CONTINUE 

WRITE ( * , 801) 

READ ( * , * ) MBASE 


DO 2 0 1 = 1, MBASE 

XDATA (I) = FINDEX (I) 
YDATA (I) « FYVAL (I) 

20 CONTINUE 


CALL UN I VAR ( MBASE, NVALS, NMISS, YBAR, YVAR, YSUM ) 


IF ( NVALS .GT. 0 ) THEN 

WRITE (OUNIT, 2001) FXVAL (1), FXVAL (MBASE) , NVALS, NMISS, 

YBAR, YVAR, YSUM 

WRITE (*,2001) FXVAL ( 1 ) , FXVAL (MBASE ), NVALS, NMISS, YBAR, 
YVAR, YSUM 

ELSE 

WRITE (*, 2000) 

RETURN 

END IF 


WRITE (*, 802) 

READ ( * , * ) ITYPE 

IF ( ITYPE .EQ. 1 ) GOTO 1 

WRITE (OUNIT, 2002) 

WRITE (*,2002) 

INITIALIZE AND COMPUTE FIXED VALUES; BEGIN PROGRESSIVE VALUE LOOP. 

BARM = YBAR 
M = FLOAT (MBASE) 

GAMOLD = FMISS 
J1 = MBASE + 1 
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J2 = NCASE 
C 

C FILL PROGRESSIVE WORK ARRAYS, COMPUTE PROGRESSIVE STATISTICS, AND 
C OUTPUT RESULTS. 

C 

C 

DO 22 JJ - Jl, J2 

XDATA(JJ) = FINDEX(JJ) 

YDATA(JJ) = FYVAL(JJ) 

CALL UNIVAR ( JJ, NVALS, NMISS, YBAR, YVAR, YSUM ) 

COMPUTE GAMMA. SEE TEXT FOR EQUATION DETAILS. 


IF ( NVALS .GT. 0 ) THEN 

BARMJ « YBAR 
J - FLOAT (JJ) - M 

T1 = (BARMJ * <M + J) ) - (BARM * M) 
T2 « ALOG (BARMJ / BARM) 

T3 = (BARMJ - BARM) * J 
Q = EXP (T1 * T2 - T3 ) 

GA.MMA = 1. / (1. + SQRT(Q)) 

D GAMMA = GAMOLD - GAMMA 


IF ( GAMOLD .EQ. FMISS ) THEN 

WRITE (OUNIT, 2003) JJ, FXVAL(JJ), YDATA(JJ), BARMJ, YVAR 

YSUM, GAMMA, CFLAG 

WRITE (*,2003) JJ, FXVAL(JJ), YDATA(JJ), BARMJ, YVAR, 
YSUM, GAMMA, CFLAG 

ELSE 

WRITE (OUNIT, 2003) JJ, FXVAL(JJ), YDATA(JJ), BARMJ, YVAR 
YSUM, GAMMA, CFLAG, D GAMMA 
WRITE (*, 2003) JJ, FXVAL(JJ), YDATA(JJ), BARMJ, YVAR, 
YSUM, GAMMA, CFLAG, D GAMMA 

END IF 


GAMOLD = GAMMA 
ELSE 

WRITE (OUNIT, 2004) JJ, FXVAL(JJ), YDATA(JJ) 
WRITE (*, 2004) JJ, FXVAL(JJ), YDATA(JJ) 
GAMOLD = FMISS 

END IF 

C 

22 CONTINUE 


+ 

+ 


+ 

+ 
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WRITE (*, 802) 

READ ( * , * ) ITYPE 

IF ( ITYPE .EQ. 1 ) GOTO 1 


801 FORMAT (/,' ENTER NUMBER OF VALUES IN BASE PERIOD:') 

802 FORMAT (/,' ENTER 0=CONTINUE, 1=CHANGE NUMBER OF BASE VALUES:') 

2000 FORMAT (' ALL VALUES IN BASE PERIOD MISSING') 

2001 FORMAT (/, 'BASE PERIOD PARAMETERS:',/, 

+ 'X-VALUE RANGE = ' , F5 . 0, ' , F5 . 0, /, 

+ 'NUMBER OF VALUES = ',13,/, 

+ 'NUMBER OF MISSING VALUES = ',13,/, 

+ 'BASE MEAN = ',F8.3,/, 

+ 'BASE VARIANCE = ',F8.3,/, 

+ 'BASE SUM = ' ,F6.0, /) 

2002 FORMAT (/, 'PROGRESSIVE PARAMETERS:',/, 

+ ' X Y 

+ ' DELTA',/, 

+ ' INDEX OBSERVATIONS MEAN VARIANCE SUM GAMMA' , 

+ ' GAMMA',/, 

+ ' ======' ) 

2003 FORMAT (15, 2X, F8.3, IX, F7 .0, 2 (2X,F8.3) , 2X, f6.0, 2X, F6 . 3, 2X, F6 . 3 ) 

2004 FORMAT (I5,2X,F8.3, IX, F7 . 0 ) 


RETURN 

END 

C 

SUBROUTINE SCHI 
C 

C SUBROUTINE TO COMPUTE CHANGE IN CHI-SQUARE DEGREES OF FREEDOM. 

C 

C VARIABLES USED IN SUBROUTINE NOT DESCRIBED IN MAIN PROGRAM: 

C 'CHANGE' PERCENT CHANGE OF ' OLDPROD' TO 'PRODZ'. 

C 'D' GAMMA FUNCTION CONSTANT. 

C ' DGAMMA' CHANGE IN GAMMA FROM PREVIOUS VALUE. 

C 'GAMMA' GAMMA VALUE. 

C 'GAMNUM' VALUE OF THE GAMMA FUNCTION. 

C ' GAMOLD' PREVIOUS VALUE OF GAMMA. 

C 'J' PROGRESSIVE VALUE INDEX. 

C 'JJ' DO LOOP COUNTER. 

C 'Jl' STARTING VALUE FOR DO LOOP COMPUTING PROGRESSIVE VALUES. 

C ' J2' ENDING VALUE FOR DO LOOP COMPUTING PROGRESSIVE VALUES. 

C 'M' REAL VALUE OF NUMBER OF BASE VALUES. 

C 'MBASE' NUMBER OF BASE PERIOD VALUES. 

C 'NUM' BASE PERIOD MEAN. 

C ' NUMJ' PROGRESSIVE MEAN. 

C 'OLDPROD' PREVIOUS VALUE OF 'PRODZ'. 

C 'PRODZ' PRODUCT OF GAMMA FUNCTION EULER RELATION. 

C 'Q' COMBINED LIKELIHOOD RATIOS. 

C ' SUMCHI2' PROGRESSIVE SUM OF CHI-SQUARE VALUES. 

C 'Tl' WORK VARIABLE; FIRST TERM IN EITHER Q EQUATION. 
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C 

C ' T2' 

C 

C ' T3' 

C ' Z' 

C 

PARAMETER (NDIM=150) 

COMMON /WORK/ FXVAL (NDIM) , FYVAL (NDIM) , XDATA (NDIM) , 

+ YDATA (NDIM) , FINDEX (NDIM) , FCHI2 (NDIM) , NCASE, 

+ OUNIT, FMISS 

REAL J, M, NUM, NUMJ 
INTEGER OUNIT 

D = .5772156649 

FILL WORK ARRAYS WITH BASE PERIOD VALUES, COMPUTE BASE PERIOD 
STATISTICS, AND OUTPUT RESULTS. 

1 CONTINUE 

WRITE (*, 801) 

READ ( * , * ) MBASE 
C 

DO 20 I = 1, MBASE 

XDATA (I) = FLOAT (I) 

YDATA (I) = FCHI2 (I) 

20 CONTINUE 

C 

CALL UNIVAR ( MBASE, NVALS, NMISS, YBAR, YVAR, YSUM ) 

C 

IF ( NVALS .GT. 0 ) THEN 

WRITE (OUNIT, 3001) FXVAL (1), FXVAL (MBASE ) , NVALS, NMISS, YBAR 
WRITE (*,3001) FXVAL (1), FXVAL (MBASE ) , NVALS, NMISS, YBAR 

ELSE 

WRITE (*, 3000) 

RETURN 

END IF 

C 

WRITE ( *, 802) 

READ ( * , * ) I TYPE 

IF ( I TYPE .EQ. 1 ) GOTO 1 

WRITE (OUNIT, 3002) 

WRITE ( * , 3002) 

C 

C COMPUTE GAMMA FUNCTION FOR THE BASE PERIOD. 

C 10,000 ITERATIONS OF THE LOOP IS ARBITRARY, BUT DOES SEEM TO ALLOW 
C FOR REASONABLE CONVERGENCE OF THE FUNCTION. 'CHANGE' IS USED TO 
C COMPUTE THE PERCENT CHANGE FROM THE PREVIOUS FUNCTION VALUE. 


FIRST TERM IN GAMMA FUNCTION EULER RELATION. 
WORK VARIABLE; SECOND TERM IN EITHER Q EQUATION. 

SECOND TERM IN GAMMA FUNCTION EULER RELATION. 
WORK VARIABLE; THIRD TERM IN EITHER Q EQUATION. 

VALUE USED IN GAMMA FUNCTION EULER RELATION. 
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C THUS, IT CAN BE USED TO EXIT FROM A LARGE ITERATION LOOP. 

C HOWEVER, WE HAVE NOT USED THIS CRITERIA CONSISTENTLY THUS FAR, BUT 
C WILL LEAVE IT BUILT INTO THE CODE FOR NOW. 

C 

NUM = YBAR 
PRODZ = 1 . 

Z = NUM / 2. 

C 

DO 26 I = 1, 10000 

T1 = 1. + Z / FLOAT (I) 

T2 = EXP (-Z / FLOAT (I)) 

PRODZ = PRODZ * T1 * T2 

C IF( I .GT. 1 ) CHANGE = ABS ( (OLDPROD - PRODZ) / OLDPROD) 

C IF ( CHANGE . LE . l.E-5 ) GOTO 27 

C OLDPROD = PRODZ 

26 CONTINUE 

C 

27 CONTINUE 

GAMNUM =1. / (Z * EXP (D * Z) * PRODZ) 

INITIALIZE AND COMPUTE FIXED VALUES; BEGIN PROGRESSIVE VALUE LOOP. 

J1 = MBASE + 1 
J2 = NCASE 
M = FLOAT (MBASE) 

GAMOLD = FMISS 
SUMCHI2 = 0. 

FILL PROGRESSIVE WORK ARRAYS, COMPUTE PROGRESSIVE STATISTICS, AND 
OUTPUT RESULTS. 


DO 22 JJ = Jl, J2 

SUMCHI2 = SUMCHI2 + ALOG (FCHI2 ( JJ) ) 

XDATA(JJ) = FLOAT (JJ) 

YDATA(JJ) = FCHI2 (JJ) 

CALL UN I VAR ( JJ, NVALS, NMISS, YBAR, YVAR, YSUM ) 


IF ( NVALS .GT. 0 ) THEN 

COMPUTE GAMMA FUNCTION FOR PROGRESSIVE VALUES. 
SEE FURTHER DESCRIPTION ABOVE. 


NUMJ = YBAR 
J = FLOAT (JJ) - M 
PRODZ = 1. 

Z = NUMJ / 2. 


ooo o noon 
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DO 28 I = 1, 10000 

T1 = 1. + Z / FLOAT (I) 

T2 = EXP (-Z / FLOAT (I)) 

PRODZ = PRODZ * T1 * T2 

IF ( I .GT. 1 ) CHANGE = ABS ( (OLDPROD - PRODZ) 
+ / OLDPROD) 

IF ( CHANGE .LE. l.E-5 ) GOTO 29 

OLDPROD = PRODZ 


28 

CONTINUE 

29 

CONTINUE 


GAMNUMJ =1. / (Z * EXP (D * Z) * PRODZ) 

COMPUTE GAMMA. SEE TEXT FOR EQUATION DETAILS. 

T1 = (J * (NUM - NUMJ) * ALOG (2 . ) ) / 2. 
T2 = J * ALOG (GAMNUM / GAMNUMJ) 

T3 = ( (NUMJ - NUM) / 2 .) * SUMCHI2 
Q = EXP (T1 + T2 + T3) 

GAMMA = 1. / (1. + SQRT(Q)) 

DGAMMA = GAMOLD - GAMMA 


IF ( GAMOLD .EQ. FMISS ) THEN 

WRITE (OUNIT, 3003) JJ, YDATA(JJ), NUMJ, GAMMA, CFLAG 
ELSE 

WRITE (OUNIT, 3003) JJ, YDATA(JJ), NUMJ, GAMMA, CFLAG, 
+ DGAMMA 


END IF 


GAMOLD = GAMMA 
ELSE 

WRITE (OUNIT, 3004) JJ, YDATA(JJ) 
WRITE (*,3004) JJ, YDATA(JJ) 
GAMOLD = FMISS 

END IF 

C 

22 CONTINUE 


WRITE ( * , 802) 

READ ( * , * ) I TYPE 

IF ( ITYPE .EQ. 1 ) GOTO 1 
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801 FORMAT ( / , ' ENTER NUMBER OF VALUES IN BASE PERIOD:') 

802 FORMAT (/,' ENTER 0=CONTINUE, 1=CHANGE NUMBER OF BASE VALUES:') 

3000 FORMAT (' ALL VALUES IN BASE PERIOD MISSING') 

3001 FORMAT (/, 'BASE PERIOD PARAMETERS:',/, 

+ 'X-VALUE RANGE = ' , F5 . 0, ' -' , F5 . 0, /, 

+ 'NUMBER OF VALUES = ',13,/, 

+ ' NUMBER OF MISSING VALUES = ',13,/, 

+ 'BASE DEGREES OF FREEDOM = ',F8.3,/) 

3002 FORMAT (/,' PROGRESSIVE PARAMETERS:',/, 

+ ' CHI-SQUARE DELTA',/, 

+ 'INDEX OBSERVATIONS MEAN GAMMA GAMMA',/, 

3003 FORMAT (I5,2X,F10.3,2X,F8.3,2X,F6.3,2X,F6.3) 

3004 FORMAT (I5,2X,F10.3) 

RETURN 

END 

C 

SUBROUTINE SLINEAR 
C 

C SUBROUTINE TO COMPUTE CHANGE IN LINEAR TREND. 

C 

C VARIABLES USED IN SUBROUTINE NOT DESCRIBED IN MAIN PROGRAM: 

C 'AM' BASE PERIOD INTERCEPT FROM LLS REGRESSION. 

C 'BM' BASE PERIOD SLOPE FROM LLS REGRESSION. 

C 'DM' CONSTANT USED IN Q (M) LIKELIHOOD RATIO. 

C ' DMJ' CONSTANT USED IN Q (M+ J) LIKELIHOOD RATIO. 

C 'DENOM' DENOMINATOR IN 'T21' AND 'T22' TERMS. 

C ' DGAMMA' CHANGE IN GAMMA FROM PREVIOUS VALUE. 

C 'GAMMA' GAMMA VALUE. 

C 'GAMOLD' PREVIOUS VALUE OF GAMMA. 

C 'I' DO LOOP COUNTER. 

C 'J' PROGRESSIVE VALUE INDEX. 

C 'JJ' DO LOOP COUNTER. 

C 'Jl' STARTING VALUE FOR DO LOOP COMPUTING PROGRESSIVE VALUES. 

C ' J2' ENDING VALUE FOR DO LOOP COMPUTING PROGRESSIVE VALUES. 

C ' LNL1' NATURAL LOG LIKELIHOOD 1. 

C ' LNL2' NATURAL LOG LIKELIHOOD 2. 

C ' LNL3' NATURAL LOG LIKELIHOOD 3. 

C ' LNL4 ' NATURAL LOG LIKELIHOOD 4. 

C 'M' REAL VALUE OF NUMBER OF BASE VALUES. 

C 'MBASE' NUMBER OF BASE PERIOD VALUES. 

C ' PREDM' PREDICTED VALUES FROM REGRESSION EQUATION USING 'M' VALUES. 

C ' PREDMJ' PREDICTED VALUES FROM REGRESSION EQUATION USING 'M+J' VALUES. 

C 'Q' RATIO OF QMJ/QM. 

C 'QM' BASE PERIOD LIKELIHOOD RATIO. 

C ' QMJ' PROGRESSIVE PERIOD LIKELIHOOD RATIO. 

C ' RESIDM' ARRAY OF BASE PERIOD RESIDUALS FROM LLS REGRESSION. 

C ' RESIDMJ' ARRAY OF PROGRESSIVE VALUE RESIDUALS FROM LLS REGRESSION. 

C ' SSM' SUM OF SQUARES FOR BASE PERIOD. 

C ' SSMJ' SUM OF SQUARES FOR 'M+J' VALUES. 
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C ' SUMT2 ' SUM OF 'T2' TERM IN ' LNL1' AND ' LNL4' FORMULAE. 

C ' SUMT2 1 ' SUM OF 'T2' TERM IN ' LNL2 ' FORMULAE. 

C ' SUMT22' SUM OF 'T2' TERM IN ' LNL3' FORMULAE. 

C 'T' REAL VALUE OF BASE PERIOD LOOP INDEX. 

C ' Tl' WORK VARIABLE FOR FIRST TERM IN ' LNL1' AND ' LNL4' FORMULAE. 

C ' T2' WORK VARIABLE FOR SECOND TERM IN ' LNL1' AND ' LNL4' FORMULAE. 

C 'Til' WORK VARIABLE FOR FIRST TERM IN ' LNL2' FORMULAE. 

C ' T12' WORK VARIABLE FOR FIRST TERM IN ' LNL3' FORMULAE. 

C ' T21' WORK VARIABLE FOR SECOND TERM IN ' LNL2 ' FORMULAE. 

C ' T22' WORK VARIABLE FOR SECOND TERM IN ' LNL3' FORMULAE. 

C 

PARAMETER (NDIM=150) 

COMMON /WORK/ FXVAL (NDIM) , FYVAL (NDIM) , XDATA (NDIM) , 

+ YDATA (NDIM) , F INDEX (NDIM) , FCHI2 (NDIM) , NCASE, 

+ OUNIT, FMISS 

REAL RESIDM (NDIM), RES IDMJ (NDIM) , PREDM (NDIM) , PREDMJ (NDIM) 

REAL M, J, LNL1 , LNL2, LNL3, LNL4 
INTEGER OUNIT 

FILL WORK ARRAYS WITH BASE PERIOD VALUES, COMPUTE BASE PERIOD 
STATISTICS, AND OUTPUT RESULTS. 

1 CONTINUE 

WRITE (*, 801) 

READ ( * , * ) MBASE 
C 

DO 20 I = 1, NCASE 

XDATA (I) = F INDEX ( I ) 

YDATA (I) = FYVAL (I) 

20 CONTINUE 

C 

CALL UNIVAR ( MBASE, NVALS, NMISS, YBAR, YVAR, YSUM ) 

C 

IF ( NVALS .GT. 0 ) THEN 

CALL RCOEFF ( MBASE, AM, BM ) 

WRITE (OUNIT, 4001) FXVAL (1), FXVAL (MBASE) , NVALS, NMISS, 

YBAR, YVAR, BM 

WRITE (*, 4001) FXVAL (1), FXVAL (MBASE) , NVALS, NMISS, YBAR, 
YVAR, BM 

ELSE 

WRITE (*, 4000) 

RETURN 

END IF 

C 

WRITE (*, 802) 

READ ( * , * ) I TYPE 

IF ( I TYPE .EQ. 1 ) GOTO 1 
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WRITE (OUNIT, 4002) 

WRITE (*, 4002) 

M = FLOAT (MBASE) 

DM = (M * (M* *2 - 1 . ) ) / 12 . 

SSM = 0. 

c 

DO 22 I = 1, MBASE 
C 

IF ( YDATA(I) .NE. FMISS ) THEN 

SSM = SSM + (YDATA(I) - YBAR) **2 
END IF 

c 

22 CONTINUE 

C 

C 

DO 23 I = 1, NCASE 
C 

IF ( YDATA(I) .NE. FMISS ) THEN 

PREDM(I) = (XDATA(I) * BM + AM) 

RESIDM(I) = (YDATA(I) -PREDM(I)) 

ELSE 

PREDM(I) = FMISS 
RESIDM (I ) = FMISS 

END IF 

C 

23 CONTINUE 

C 

SUMT2 = 0. 

T1 = (-1.) * (M - 2.) / (2. * (SSM - (BM* *2 * DM))) 

C 

DO 24 I = 1, MBASE 
C 

IF ( YDATA(I) .NE. FMISS ) THEN 
T = FLOAT (I) 

T2 = RESIDM(I) **2 / ( ( (M + 1 . ) / M) 

+ + (T - ( (M + 1.) / 2.) )**2 / DM) 

SUMT2 = SUMT2 + T2 

END IF 

C 

24 CONTINUE 

LNL1 = T1 * SUMT2 


INITIALIZE AND COMPUTE FIXED VALUES; BEGIN PROGRESSIVE VALUE LOOP. 
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C 

GAMOLD = FMISS 
J1 = MBASE + 1 
J2 = NCASE 
C 

C FILL PROGRESSIVE WORK ARRAYS, COMPUTE PROGRESSIVE STATISTICS, AND 
C OUTPUT RESULTS . 

C 

C 

DO 26 JJ = Jl, J2 

XDATA(JJ) = FINDEX(JJ) 

YDATA(JJ) = FYVAL(JJ) 

CALL UNIVARt JJ, NVALS, NMISS, YBAR, YVAR, YSUM ) 

IF ( YDATA(JJ) .EQ. FMISS ) NVALS = 0 
COMPUTE EQUATION TERMS AND GAMMA. SEE TEXT FOR EQUATION DETAILS. 


IF ( NVALS -GT. 0 ) THEN 

CALL RCOEFF ( JJ, AMJ, BMJ ) 


SSMJ = 0. 

J = FLOAT (JJ) - M 

DMJ = ( (M + J) * ( (M + J) * *2 - 1.)) / 12. 



DO 30 I = 1, JJ 


SSMJ = SSMJ + (YDATA(I) - YBAR)**2 
PREDMJ ( I ) = (XDATA(I) * BMJ + AMJ) 
RESIDMJ ( I ) = (YDATA(I) - PREDMJ(I)) 

30 

CONTINUE 


SUMT2 * 0 . 

T1 = (-1.) * (M - 2.) / (2. * (SSM - (BMJ* *2 


DO 32 I = 1, MBASE 


T = FLOAT (I) 

+ 

T2 = RESIDMJ (I ) **2 / ( ( (M + 1) / M) + 

(T - <(M + 1.) / 2.))**2 / DM) 


SUMT2 = SUMT2 + T2 

32 

C 

CONTINUE 


LNL4 = T1 * SUMT2 


SUMT21 = 0. 
SUMT22 = 0. 
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C 



Til = 

(-1.) 

* (M + J - 

2.) 

+ 


/ (2. 

* (SSMJ - 

(BMJ* *2 


T12 = 

(-1.) 

* (M + J - 

2.) 

+ 


/ (2. 

* (SSMJ - 

(BM* *2 * 


DO 34 

1 = 1, 

JJ 



* DMJ) ) ) 
DMJ) ) ) 


T = FLOAT (I) 


DENOM = ( (M + J + 1 . ) / (M + J) ) 

+ + ( (T — ( (M + J + 1.) / 2.))**2 / DMJ) 

T21 = RESIDMJ (I) **2 / DENOM 
122 = RESIDM (I ) **2 / DENOM 

SUMT21 = SUMT21 + T21 

SUMT22 = SUMT22 + T22 

34 CONTINUE 

C 

LNL2 = Til * SUMT21 
LNL3 = T12 * SUMT22 


QM = EXP (LNL4 - LNL1) 

QMJ = EXP (LNL2 - LNL3) 

Q = QMJ / QM 

GAMMA = 1. / (1. + SQRT(Q)) 
D GAMMA = GAMOLD - GAMMA 


IF ( GAMOLD .EQ. FMISS ) THEN 

WRITE (OUNIT, 4003) JJ, FXVAL(JJ), YDATA(JJ), YBAR, YVAR 
+ BMJ, GAMMA 

WRITE (*, 4003) JJ, FXVAL(JJ), YDATA(JJ), YBAR, YVAR, 

+ BMJ, GAMMA 


ELSE 


+ 

+ 


WRITE (OUNIT, 4003) JJ, FXVAL(JJ), YDATA(JJ), YBAR, YVAR 

BMJ, GAMMA, D GAMMA 

WRITE (*, 4003) JJ, FXVAL(JJ), YDATA(JJ), YBAR, YVAR, 

BMJ, GAMMA, DGAMMA 


C 


END IF 


GAMOLD = GAMMA 


ELSE 


WRITE (OUNIT, 4004) JJ, FXVAL(JJ), YDATA(JJ) 
WRITE (*, 4004) JJ, FXVAL(JJ), YDATA(JJ) 

END IF 

C 

26 CONTINUE 
C 
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WRITE (*, 802) 

READ ( * , * ) ITYPE 

IF ( ITYPE .EQ. 1 ) GOTO 1 

801 FORMAT (/,' ENTER NUMBER OF VALUES IN BASE PERIOD:') 

802 FORMAT (/,' ENTER 0=CONTINUE, 1=CHANGE NUMBER OF BASE VALUES:') 

4000 FORMAT (' ALL VALUES IN BASE PERIOD MISSING') 

4001 FORMAT (/, 'BASE PERIOD PARAMETERS:',/, 

+ 'X-VALUE RANGE = ' , F5 . 0, ' -' , F5 . 0, /, 

+ 'NUMBER OF VALUES = ',13,/, 

+ 'NUMBER OF MISSING VALUES = ',13,/, 

+ 'BASE MEAN = ',F8.3,/, 

+ 'BASE VARIANCE = ',F8.3,/, 

+ 'BASE SLOPE = ' , F6 . 3, /) 

4002 FORMAT (/,' PROGRESSIVE PARAMETERS:',/, 

+ ' X Y 

+ ' DELTA' , /, 

+ ' INDEX OBSERVATIONS MEAN VARIANCE 

+ ' GAMMA',/, 

+ ' ======' ) 

4003 FORMAT (15, 2X,F8 .3, IX, F7 .2, 2 (2X, F8 . 3) , 2 (2X,F6.3) , 2X,F6.3) 

4004 FORMAT (15, 2X,F8 . 3, IX, F7 .2) 


f 

/ 

SLOPE GAMMA' , 



C- 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


RETURN 

END 


SUBROUTINE RCOEFF ( N, B0, B1 ) 

SUBROUTINE TO COMPUTE LINEAR LEAST SQUARES (LLS) INTERCEPT AND SLOPE. 

VARIABLES USED IN SUBROUTINE NOT DESCRIBED IN MAIN PROGRAM: 

'I' DO LOOP COUNTER. 

'N' NUMBER OF VALUES. 

'NN' COUNTER FOR NON-MISSING VALUES. 

' SX2' SUMS OF SQUARED DEVIATIONS FROM THE MEAN X-VALUE. 

' SXY' SUM OF THE CROSS PRODUCTS OF DEVIATIONS. 

' SUMX' SUM OF NON-MISSING X-VALUES. 

'SUMY' SUM OF NON-MISSING Y-VALUES . 

' SUMSX2' SUM OF X-VALUES SQUARED. 

' SUMSXY' SUM OF X TIMES Y SQUARED. 

' XBAR' MEAN OF NON-MISSING X-VALUES. 

' XBAR2 ' NUMBER OF NON-MISSING VALUES TIMES 'XBAR' SQUARED. 

' XWORK' WORK ARRAY FOR X-VALUES. 

' XYN' 'XBAR' TIMES ' YBAR' TIMES NUMBER OF NON-MISSING VALUES. 
'YBAR' MEAN OF NON-MISSING Y-VALUES. 

' YWORK' WORK ARRAY FOR Y-VALUES. 

PARAMETER (NDIM=150) 

COMMON /WORK/ FXVAL (NDIM) , FYVAL (NDIM) , XDATA(NDIM), 

+ YDATA (NDIM) , F INDEX (NDIM ) , FCHI2 (NDIM) , NCASE , 
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+ OUNIT, FMISS 

REAL XWORK (NDIM) , YWORK (NDIM) 

INTEGER OUNIT 

SUMX - 0. 

SUMY - 0. 

NN = 0 
C 

C FILL WORK ARRAYS WITH NON-MISSING VALUES. 
C 

c 

DO 15 I = 1, N 


IF ( YDATA(I) .GT. FMISS ) THEN 

NN = NN + 1 

XWORK (NN) = XDATA(I) 

YWORK (NN) = YDATA(I) 

END IF 

C 

15 CONTINUE 

C 

C 

C COMPUTE SUMS OF X- AND Y-VALUES . 

C 

C 

DO 20 I = 1, NN 

SUMX = SUMX + XWORK (I) 

SUMY = SUMY + YWORK (I) 

20 CONTINUE 

C 

C 

C COMPUTE MEANS OF X' AND Y-VALUES. 

C 

XBAR = SUMX / FLOAT (NN) 

YBAR = SUMY / FLOAT (NN) 

SUMSX2 = 0. 

SUMSXY = 0. 

C 

C COMPUTE SUMS OF SQUARED DEVIATIONS AND CROSS PRODUCTS . 
C 

C 

DO 22 I = 1, NN 

SUMSX2 = SUMSX2 + XWORK(I)**2 

SUMSXY = SUMSXY 4 - (XWORK (I) * YWORK (I)) 

22 CONTINUE 


XBAR2 = XBAR* *2 * FLOAT (NN) 
XYN = XBAR * YBAR * FLOAT (NN) 
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SX2 = SUMSX2 - XBAR2 
SXY = SUMSXY - XYN 

COMPUTE COEFFICIENTS. 

B1 = SXY / SX2 

BO = YBAR - B1 * XBAR 

RETURN 

END 


SUBROUTINE UNIVAR( N, NVALS, NMISS, YBAR, YVAR, YSUM ) 

SUBROUTINE TO COMPUTE LINEAR LEAST SQUARES (LLS) INTERCEPT AND SLOPE. 

VARIABLES USED IN SUBROUTINE NOT DESCRIBED IN MAIN PROGRAM: 

'I' DO LOOP COUNTER. 

'N' NUMBER OF VALUES. 

' SUMS2' SUM OF SQUARED DIFFERENCES. 

' YWORK' WORK ARRAY FOR NON-MISSING Y-VALUES. 

PARAMETER (NDIM=150) 

COMMON /WORK/ FXVAL(NDIM), FYVAL(NDIM), XDATA(NDIM), 

+ YDATA (NDIM) , FINDEX (NDIM) , FCHI2 (NDIM) , NCASE, 

+ OUNIT, FMISS 

REAL YWORK (NDIM) 

INTEGER OUNIT 

YSUM = 0. 

NVALS = 0 
NMISS = 0 

FILL WORK ARRAY WITH NON-MISSING VALUES. 


DO 15 I = 1, N 

IF ( YDATA (I) .GT. FMISS ) THEN 

NVALS = NVALS + 1 
YWORK (NVALS) = YDATA (I) 

ELSE 

NMISS = NMISS + 1 
END IF 
CONTINUE 


COMPUTE SUM OF Y-VALUES. 


DO 20 I = 1, NVALS 
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YSUM = YSUM + YWORK(I) 

20 CONTINUE 

COMPUTE MEAN OF Y-VALUES . 

YBAR = YSUM / FLOAT (NVALS) 

COMPUTE SUM OF SQUARED DIFFERENCES. 

SUMS2 = 0. 

DO 22 I = 1, NVALS 

SUMS 2 = SUMS 2 + (YWORK(I) - YBAR)* *2 
22 CONTINUE 

COMPUTE VARIANCE OF Y-VALUES. 

YVAR = SUMS 2 / (FLOAT (NVALS) - 1) 

RETURN 

END 

C 


